Regressions for T0 Stress

Set Gobal Chunk Options & Load libraries

• WARNING: Setting cache = TRUE can cause random errors in the code (like variables that cannot be found suddenly). If error arises, set cache=FALSE and rerun.

knitr::opts_chunk$set(echo = T,warning = FALSE, cache = FALSE, collapse = T)
getwd()
#library(nlme) #mixed modeling - delete if not used
library(tidyverse)
library(plyr)#revalue() factors and reorder them! dplyr::recode not working
library(dplyr)  #summarise_at()
library(ggplot2) #graphing
library(ggfortify) #plot multiple plots from model output into one panel!
#library(lattice) #delete if not used
#library(gtable)
library(grid) #tables?
library(gridExtra)#masks combine() from dplyr
library(pander) #table output
library(Publish) #delete if not used
library(ggpubr) #combine plots
library(kableExtra)

## Warning in !is.null(rmarkdown::metadata$output) && rmarkdown::metadata$output
## %in% : 'length(x) = 2 > 1' in coercion to 'logical(1)'

library(devtools) #
library(broom) #
library(purrr)#

#Modeling libraries
library(car)#sophisticated anova() modeling
library(rstatix)#anova_test() - masks mutate() and desc() from plyr package!
library(sjstats)#anova_stats()
library(pwr)#unknown function but the sjstats package must use it to calculate Cohen's D!
#library(MDmisc)#download from GitHub - 
library(ggiraphExtra)
library(ggcorrplot)# modern correlation matrix plot
library(jtools)# regression tables etc >summ()

Import Data & Set Factor Levels

• Question: How can you convert a bunch of columns to factors simultaneously?

#rm(Recovery)
D<-read.csv("T0_Stress_R.csv")

# Set factors (no longer automatic)
D$Moss<-as.factor(D$Moss)
D$Microsite<-as.factor(D$Microsite)

# Order factor levels
D$Habitat<-factor(D$Habitat, levels= c("Interspace", "Dripline", "Canopy"))
D$Topography<-factor(D$Topography, levels=c("S", "F", "N"))
D$Veg.Type<-factor(D$Veg.Type, levels=c("CB", "BB", "PJ"))
D$Veg.Type<-plyr::revalue(D$Veg.Type, c("CB" = "Low-scrubland", "BB" = "Mid-shrubland", "PJ" ="High-woodland"))


D$Site<- factor(D$Site, levels=c("1","2","3"))
#levels(D$Site)<- list("1" = "Low-scrubland", "2"= "Mid-shrubland", "3" = "High-woodland")

#Remove other species (Syntrichia ruralis) & drop microsites with PAR sensors, not used for tissue sampling (35 & 53).
D<-droplevels(D[!c(D$Microsite=="PJ-N-91" | D$Microsite =="PJ-S-67a" | 
                    D$Microsite=="BB-N-35" | D$Microsite=="BB-S-53" |D$Microsite=="CB-F-1b"| D$Microsite =="CB-F-5b"),])#35 & 53 PAR sensor plots have no stress data but DO have habitat data.


length(levels(D$Microsite))#Down to 94 units
## [1] 92

MUTATE VARIABLES

Should I model with shade as a percent? I know percentages can have strange properties…but it is still a ratio by definition with upper and lower limits.

#Convert Shade.Index to Percent
#D<-D %>% mutate(Shade.Index = (Shade.Index/28)*100)
#Convert PDIR to units in table (kJ/cm²/yr), was in joules
#D<-D %>% mutate(PDIR= PDIR*1000)
D<- D %>% mutate(Shade.Time=(Shade.Index/28)*100)
summary(D)
##       Time        Moss     Moss.number      Moss.sub            Microsite 
##  Min.   :0   10     : 1   Min.   : 1.00   Length:92          BB-F-37 : 1  
##  1st Qu.:0   11     : 1   1st Qu.:23.75   Class :character   BB-F-38a: 1  
##  Median :0   12     : 1   Median :48.50   Mode  :character   BB-F-39a: 1  
##  Mean   :0   13a    : 1   Mean   :48.08                      BB-F-40a: 1  
##  3rd Qu.:0   14a    : 1   3rd Qu.:72.25                      BB-F-41a: 1  
##  Max.   :0   15a    : 1   Max.   :96.00                      BB-F-42a: 1  
##              (Other):86                                      (Other) :86  
##  Old.Habitat              Habitat       Index         Elevation    Site  
##  Length:92          Interspace:17   Min.   : 3.00   Min.   : 890   1:24  
##  Class :character   Dripline  :27   1st Qu.: 8.00   1st Qu.: 890   2:34  
##  Mode  :character   Canopy    :48   Median :14.00   Median :1670   3:34  
##                                     Mean   :14.67   Mean   :1614         
##                                     3rd Qu.:20.50   3rd Qu.:2070         
##                                     Max.   :26.00   Max.   :2070         
##                                                                          
##           Veg.Type  Topography  Shade.Index    Summer.Shade.Index
##  Low-scrubland:24   S:23       Min.   : 6.00   Min.   : 2.00     
##  Mid-shrubland:34   F:35       1st Qu.:12.00   1st Qu.: 9.00     
##  High-woodland:34   N:34       Median :20.00   Median :23.00     
##                                Mean   :18.09   Mean   :21.21     
##                                3rd Qu.:23.00   3rd Qu.:32.25     
##                                Max.   :27.00   Max.   :42.00     
##                                                                  
##       PDIR            Slope            Aspect        VegPercent    
##  Min.   :0.7717   Min.   : 0.000   Min.   :  0.0   Min.   :  2.78  
##  1st Qu.:0.9036   1st Qu.: 2.000   1st Qu.: 90.0   1st Qu.: 24.30  
##  Median :0.9846   Median : 7.000   Median :200.0   Median : 38.89  
##  Mean   :0.9478   Mean   : 7.587   Mean   :193.2   Mean   : 41.46  
##  3rd Qu.:0.9997   3rd Qu.:12.000   3rd Qu.:270.0   3rd Qu.: 61.11  
##  Max.   :1.0496   Max.   :18.000   Max.   :357.0   Max.   :100.00  
##                                    NA's   :4                       
##     SolarVeg          Camera           Pair       Notes.Urgent      
##  Min.   :  4.17   Min.   : 1.00   Min.   : 1.00   Length:92         
##  1st Qu.: 29.17   1st Qu.: 6.25   1st Qu.:10.25   Class :character  
##  Median : 45.83   Median :11.50   Median :19.50   Mode  :character  
##  Mean   : 46.11   Mean   :11.73   Mean   :19.37                     
##  3rd Qu.: 63.54   3rd Qu.:16.75   3rd Qu.:28.75                     
##  Max.   :100.00   Max.   :23.00   Max.   :37.00                     
##                   NA's   :70      NA's   :54                        
##   Assay.Set             Assay            Clip             Fo        
##  Length:92          Min.   :1.000   Min.   : 1.00   Min.   : 30.00  
##  Class :character   1st Qu.:1.000   1st Qu.: 6.75   1st Qu.: 79.00  
##  Mode  :character   Median :2.000   Median :12.50   Median : 95.50  
##                     Mean   :2.467   Mean   :12.42   Mean   : 98.72  
##                     3rd Qu.:3.000   3rd Qu.:18.25   3rd Qu.:119.00  
##                     Max.   :4.000   Max.   :24.00   Max.   :174.00  
##                                                                     
##        Fm              Fv             Fv.Fm              Fs        
##  Min.   : 52.0   Min.   : 12.00   Min.   :0.1300   Min.   : 29.00  
##  1st Qu.:128.8   1st Qu.: 44.75   1st Qu.:0.3490   1st Qu.: 75.25  
##  Median :182.5   Median : 91.00   Median :0.4770   Median : 91.50  
##  Mean   :193.9   Mean   : 95.16   Mean   :0.4529   Mean   : 95.90  
##  3rd Qu.:248.2   3rd Qu.:132.00   3rd Qu.:0.5673   3rd Qu.:115.00  
##  Max.   :419.0   Max.   :265.00   Max.   :0.7370   Max.   :190.00  
##                                                                    
##       Fm.L            Fo.L             Fv.L           Fv.Fm.L      
##  Min.   : 46.0   Min.   : 17.00   Min.   : 15.00   Min.   :0.2230  
##  1st Qu.:100.8   1st Qu.: 51.75   1st Qu.: 42.00   1st Qu.:0.4047  
##  Median :136.5   Median : 70.00   Median : 61.00   Median :0.4890  
##  Mean   :141.1   Mean   : 70.58   Mean   : 70.49   Mean   :0.4877  
##  3rd Qu.:174.5   3rd Qu.: 90.00   3rd Qu.: 95.25   3rd Qu.:0.5790  
##  Max.   :319.0   Max.   :149.00   Max.   :211.00   Max.   :0.7360  
##                                                                    
##     PhiPSII             qP              qNP               NPQ          
##  Min.   :0.0430   Min.   :0.1330   Min.   :-0.8570   Min.   :-0.19100  
##  1st Qu.:0.1940   1st Qu.:0.4968   1st Qu.: 0.1580   1st Qu.: 0.08375  
##  Median :0.2685   Median :0.5820   Median : 0.4100   Median : 0.31850  
##  Mean   :0.2943   Mean   :0.5847   Mean   : 0.3634   Mean   : 0.43336  
##  3rd Qu.:0.3680   3rd Qu.:0.6783   3rd Qu.: 0.5883   3rd Qu.: 0.69700  
##  Max.   :0.6930   Max.   :0.9460   Max.   : 0.7970   Max.   : 2.54200  
##                                                                        
##    Shrubindex       Topoindex         Siteindex        Multibuffer     
##  Min.   :-1.000   Min.   :-1.0000   Min.   :-1.0000   Min.   :-2.0000  
##  1st Qu.: 0.000   1st Qu.:-0.2500   1st Qu.:-1.0000   1st Qu.: 0.0000  
##  Median : 1.000   Median : 0.0000   Median : 0.0000   Median : 1.0000  
##  Mean   : 0.337   Mean   : 0.1196   Mean   : 0.1087   Mean   : 0.5652  
##  3rd Qu.: 1.000   3rd Qu.: 1.0000   3rd Qu.: 1.0000   3rd Qu.: 1.0000  
##  Max.   : 1.000   Max.   : 1.0000   Max.   : 1.0000   Max.   : 3.0000  
##                                                                        
##    Shade.Time   
##  Min.   :21.43  
##  1st Qu.:42.86  
##  Median :71.43  
##  Mean   :64.60  
##  3rd Qu.:82.14  
##  Max.   :96.43  
## 

PLOT 1: Bubble Buffers Stress vs. Shade Plot

ggplot(D)+
 stat_smooth(aes(x=Shade.Index, y=Fv.Fm), formula = y ~ x, method = "lm", se = TRUE, color="grey", linetype= "solid", size=0.5)+
 geom_point(aes(x=Shade.Index, y=Fv.Fm, fill=Veg.Type, size =PDIR), shape= 21, color = "black") +
 scale_size(name = "PDIR", range = c(0.25, 6)) + 
 stat_smooth(aes(x=Shade.Index, y=Fv.Fm, color=Veg.Type), 
             formula = y ~ x, method = "lm", se = FALSE)+
 
 labs(y='\u2190 Stress (Fv/Fm)', x="Microscale shade time (% of the year)")+
 scale_y_continuous(breaks = seq(0.1, 0.75, 0.1), limits=c(0.1,0.75))+
 #annotate(geom="text",label= "R-sq = 0.082\nP < 0.0051", x=5, y=0.7)+
 #scale_color_manual(values =c("#E6AB02", "firebrick3", "dodgerblue3"))+
 scale_fill_manual(values =c("#E6AB02", "firebrick3", "dodgerblue3"), name= "Elevation\n(life zone site)", guide = guide_legend(reverse=TRUE))+ 
  scale_color_manual(values =c("#E6AB02", "firebrick3", "dodgerblue3"), name= "Life zone", guide="none")+ 
 #guides(color = "none")#guide="none")+
 theme_linedraw()+
 theme(axis.text = element_text(size=12), axis.title = element_text(size=15)) #text = element_text(size=15))#must place after chart theme

PLOT: Bubble Buffers Stress vs Shade by TOPOZONE

Not a lot of clustering (spatial autocorrelation) by topography zone, so I think I’m clear to model without blocking by Topography. * Is there a test for spatial autocorrelation?

ggplot(D)+
 geom_point(aes(x=((Shade.Index/28)*100), y=Fv.Fm, fill=interaction(Veg.Type,Topography), size = (PDIR*1000)), shape= 21, color = "black") +
 scale_size(name = "PDIR", range = c(0.05, 10)) + 
 stat_smooth(aes(x=((Shade.Index/28)*100), y=Fv.Fm, color= 
           interaction(Veg.Type,Topography)),formula = y ~ x, 
           method = "lm", se = FALSE)+
 #annotate(geom="text",label= "R-sq = 0.082\nP < 0.0051", x=5, y=0.7)+
 #scale_color_manual(values =c("#E6AB02", "firebrick3", "dodgerblue3"))
 theme_classic()

PDIR vs Stress

ggplot(D)+
 stat_smooth(formula= y~x, method="lm",
             aes(x=PDIR, y=Fv.Fm),color="grey", se=TRUE)+
 geom_point(aes(x=PDIR, y=Fv.Fm, fill=interaction(Veg.Type,Topography), 
                size = Shade.Index), shape= 21, color = "black") +
 scale_size(name = "Percent shade", range = c(1, 8)) + 
 #stat_smooth(formula = y ~ x, method = "lm", se = FALSE,
  # aes(x=PDIR, y=Fv.Fm, color= interaction(Veg.Type,Topography)))+
#geom_point(position="mean")+
 theme_classic()

#———————————

Pub 1 Analysis:

CONTINUOUS BUFFERS REGRESSION MODEL

• Microscale: sunseeker annual shade index (avoid the %time??)
• Mesoscale: PDIR (technically measured on all microsites but incorporates a larger spatial area around the moss microsite)
• Macroscale: elevation

Muliticollinearity

• Elevation has a 0.7 correlation with shade! Best to remove elevation (less precise variable) from the regression to prevent issues of multicollinearity, which inflate standard errors and can create biased predictions and unstable estimates/p-values!
• See Clark Statistics Bible or https://medium.com/analytics-vidhya/a-step-by-step-guide-for-running-a-complete-multiple-linear-regression-analysis-in-r-c08be169fe01

# Compute correlation at 2 decimal places
Buffers<-subset(D, select = c(Elevation, PDIR, Shade.Index))
Corr_matrix = round(cor(Buffers), 2)


# Compute and show the  result
ggcorrplot(Corr_matrix, hc.order = TRUE, type = "lower",
          lab = TRUE)

FULL BUFFER MODEL

• R2= 0.49 (adjusted) - .59 (unadjusted)!!
• Why does ANOVA model have way more significant factors than the linear regression model? Because the multicollinearity inflates SE of coefficients in the regression, so their p-values go way up!

Model.full<-lm(Fv.Fm ~ Elevation*PDIR*Shade.Index, data = D)
summary.lm(Model.full)#no terms are significant, but model is!
## 
## Call:
## lm(formula = Fv.Fm ~ Elevation * PDIR * Shade.Index, data = D)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.35422 -0.04914  0.00438  0.05322  0.22113 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)
## (Intercept)                 2.954e-01  1.685e+00   0.175    0.861
## Elevation                   1.914e-03  1.415e-03   1.353    0.180
## PDIR                        2.770e-01  1.770e+00   0.156    0.876
## Shade.Index                -2.541e-02  9.569e-02  -0.266    0.791
## Elevation:PDIR             -2.316e-03  1.474e-03  -1.571    0.120
## Elevation:Shade.Index      -5.336e-05  6.671e-05  -0.800    0.426
## PDIR:Shade.Index            2.625e-02  1.029e-01   0.255    0.799
## Elevation:PDIR:Shade.Index  6.893e-05  7.060e-05   0.976    0.332
## 
## Residual standard error: 0.1041 on 84 degrees of freedom
## Multiple R-squared:  0.5546, Adjusted R-squared:  0.5175 
## F-statistic: 14.94 on 7 and 84 DF,  p-value: 1.615e-12

#ANOVA of continuous variables show significance better:
summary.aov(Model.full)#Elev-shade interaction ***!!!!
##                            Df Sum Sq Mean Sq F value   Pr(>F)    
## Elevation                   1 0.0402  0.0402   3.710  0.05747 .  
## PDIR                        1 0.4068  0.4068  37.514 2.81e-08 ***
## Shade.Index                 1 0.4585  0.4585  42.288 5.36e-09 ***
## Elevation:PDIR              1 0.0034  0.0034   0.314  0.57687    
## Elevation:Shade.Index       1 0.1193  0.1193  11.006  0.00134 ** 
## PDIR:Shade.Index            1 0.0956  0.0956   8.813  0.00390 ** 
## Elevation:PDIR:Shade.Index  1 0.0103  0.0103   0.953  0.33172    
## Residuals                  84 0.9108  0.0108                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Variation Inflation Factor Test (VIF)
vif(Model.full, type = "predictor")# remove interactions before testing VIF!!
## GVIFs computed for predictors
##             GVIF Df GVIF^(1/(2*Df))         Interacts With Other Predictors
## Elevation      1  7               1      PDIR, Shade.Index             --  
## PDIR           1  7               1 Elevation, Shade.Index             --  
## Shade.Index    1  7               1        Elevation, PDIR             --
#what do these results mean??

Maine Effects (Addition) Model

#Fit Main Effects model
Model.add<-lm(Fv.Fm*100 ~ Shade.Index+PDIR+Elevation, data= D) #+PDIR+Shade.Index, data = D)
summary.lm(Model.add)#Elev-shade interaction ***!!!!
## 
## Call:
## lm(formula = Fv.Fm * 100 ~ Shade.Index + PDIR + Elevation, data = D)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -29.6381  -6.5848  -0.7143   9.1912  22.1794 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  80.930830  15.249533   5.307 8.24e-07 ***
## Shade.Index   1.730765   0.290845   5.951 5.32e-08 ***
## PDIR        -56.887801  16.155274  -3.521 0.000683 ***
## Elevation    -0.008072   0.003792  -2.129 0.036057 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.38 on 88 degrees of freedom
## Multiple R-squared:  0.4428, Adjusted R-squared:  0.4238 
## F-statistic: 23.31 on 3 and 88 DF,  p-value: 3.398e-11
#Almost X as much variation explained when Site is added to the Shade model!!!
summary.aov(Model.add)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Shade.Index  1   5626    5626  43.449 3.10e-09 ***
## PDIR         1   2843    2843  21.953 1.01e-05 ***
## Elevation    1    587     587   4.532   0.0361 *  
## Residuals   88  11395     129                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Variation Inflation Factor Test (VIF)
vif(Model.add, type = "predictor")# remove interactions before testing VIF!!
## VIFs computed for predictors
## [1] 2.067875 1.172560 2.193498
#[1] 2.281592 1.168573 2.177651 are not bad!! < 5 cutoff, elevation is highest; when n = 94
#[1] 2.067875 1.172560 2.193498 [when n = 92]

Full model no interactions

“Centering the variables is also known as standardizing the variables by subtracting the mean. This process involves calculating the mean for each continuous independent variable and then subtracting the mean from all observed values of that variable. Then, use these centered variables in your model. Most statistical software provides the feature of fitting your model using standardized variables.” - Statistics by Jim

Model.Add<-lm(Fv.Fm*100 ~ Elevation+PDIR+Shade.Index, data = D)
summary.lm(Model.Add)#Elev-shade interaction ***!!!!
## 
## Call:
## lm(formula = Fv.Fm * 100 ~ Elevation + PDIR + Shade.Index, data = D)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -29.6381  -6.5848  -0.7143   9.1912  22.1794 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  80.930830  15.249533   5.307 8.24e-07 ***
## Elevation    -0.008072   0.003792  -2.129 0.036057 *  
## PDIR        -56.887801  16.155274  -3.521 0.000683 ***
## Shade.Index   1.730765   0.290845   5.951 5.32e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.38 on 88 degrees of freedom
## Multiple R-squared:  0.4428, Adjusted R-squared:  0.4238 
## F-statistic: 23.31 on 3 and 88 DF,  p-value: 3.398e-11
#Almost X as much variation explained when Site is added to the Shade model!!!
summary.aov(Model.Add)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Elevation    1    402     402   3.107   0.0814 .  
## PDIR         1   4068    4068  31.415 2.36e-07 ***
## Shade.Index  1   4585    4585  35.412 5.32e-08 ***
## Residuals   88  11395     129                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Variation Inflation Factor Test (VIF)
Buffers.3 <- lm(Fv.Fm ~ Elevation+PDIR+Shade.Index, data= D)#how to center?
vif(Buffers.3, type = "predictor")# remove interactions before testing VIF!!
## VIFs computed for predictors
## [1] 2.193498 1.172560 2.067875
#[1] 2.281592 1.168573 2.177651 are not bad!! < 5 cutoff, elevatin is highest

Reduced Model - BEST SO FAR?

• Run reduced model with only significant interaction terms (but this changes the SDerrors and significance…suggesting multicollinearity is present.)
  • Elevation:Shade.Index
  • PDIR:Shade.Index

Model2<-lm(Fv.Fm ~ PDIR + Shade.Index, data = D)
summary.lm(Model2)# R2 = Multiple R-sq:0.4141,  Adjusted R-sq:0.4009 
## 
## Call:
## lm(formula = Fv.Fm ~ PDIR + Shade.Index, data = D)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.29752 -0.05760  0.00336  0.08733  0.22647 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.883549   0.151372   5.837 8.50e-08 ***
## PDIR        -0.699878   0.152310  -4.595 1.42e-05 ***
## Shade.Index  0.012863   0.002065   6.230 1.51e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.116 on 89 degrees of freedom
## Multiple R-squared:  0.4141, Adjusted R-squared:  0.4009 
## F-statistic: 31.45 on 2 and 89 DF,  p-value: 4.655e-11
summary.aov(Model2)# R2 = 
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## PDIR         1 0.3244  0.3244   24.10 4.12e-06 ***
## Shade.Index  1 0.5225  0.5225   38.81 1.51e-08 ***
## Residuals   89 1.1981  0.0135                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
vif(Model2)#PDIR & Shade.Index = 1.002445    1.002445 
##        PDIR Shade.Index 
##    1.002445    1.002445

PDIR * Shade Model

Because Elevation is highly correlated with shade (cor = 0.70), removing the less descriptive variable makes sense especially because elevation-life zone ANOVA was not statistically significant for explaining Fv/Fm summer stress. * The interaction of PDIR * Shade was not significant, so this term is left out.

Model.final<-lm(Fv.Fm ~ PDIR * Shade.Index, data = D)
summary.lm(Model.final)# R2 = 0.3649, P < 0.0001 DF(2,91)
## 
## Call:
## lm(formula = Fv.Fm ~ PDIR * Shade.Index, data = D)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.31212 -0.05584  0.01289  0.06602  0.22560 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       1.97825    0.52206   3.789 0.000276 ***
## PDIR             -1.86056    0.55130  -3.375 0.001100 ** 
## Shade.Index      -0.04157    0.02497  -1.665 0.099536 .  
## PDIR:Shade.Index  0.05783    0.02644   2.187 0.031401 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1136 on 88 degrees of freedom
## Multiple R-squared:  0.4443, Adjusted R-squared:  0.4254 
## F-statistic: 23.45 on 3 and 88 DF,  p-value: 3.022e-11
summary.aov(Model.final)# both are highly significant! If the variables are centered, will Shade become signif on its own?
##                  Df Sum Sq Mean Sq F value   Pr(>F)    
## PDIR              1 0.3244  0.3244  25.120 2.76e-06 ***
## Shade.Index       1 0.5225  0.5225  40.459 8.73e-09 ***
## PDIR:Shade.Index  1 0.0618  0.0618   4.783   0.0314 *  
## Residuals        88 1.1364  0.0129                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Jtools Regression Tutorial with summ() magic

I am multiplying the response stress by 100 for the tutorial! consider reporting it this way in the paper too! = Percent efficiency!!

Raw Model

library(jtools)
Model.Final<-lm(Fv.Fm ~ PDIR * Shade.Index, data = D)
summ(Model.Final, digits = 4)

Observations	92
Dependent variable	Fv.Fm
Type	OLS linear regression

F(3,88)	23.4537
R²	0.4443
Adj. R²	0.4254

	Est.	S.E.	t val.	p
(Intercept)	1.9783	0.5221	3.7893	0.0003
PDIR	-1.8606	0.5513	-3.3748	0.0011
Shade.Index	-0.0416	0.0250	-1.6647	0.0995
PDIR:Shade.Index	0.0578	0.0264	2.1869	0.0314
Standard errors: OLS

Centered & Fv/Fm Percent*100

Yes! Shade became significant as a main effect after centering (because centering removes the structural multicollinearity from the interaction term.)

summ(Model.Final, center = TRUE, confint=TRUE, digits = 4)#R2-adj = 

Observations	92
Dependent variable	Fv.Fm
Type	OLS linear regression

F(3,88)	23.4537
R²	0.4443
Adj. R²	0.4254

	Est.	2.5%	97.5%	t val.	p
(Intercept)	0.4542	0.4306	0.4778	38.2858	0.0000
PDIR	-0.8147	-1.1289	-0.5004	-5.1516	0.0000
Shade.Index	0.0132	0.0092	0.0173	6.5216	0.0000
PDIR:Shade.Index	0.0578	0.0053	0.1104	2.1869	0.0314
Standard errors: OLS; Continuous predictors are mean-centered.

#kable(Out) not working until RMarkdown runs!

_____________________

Residual Diagnostics

ggplot2:autoplot() plots all residual plots in the same window!

plot(Model.Final)

#autoplot(Model12, which = 1:4, ncol = 2, label.size = 3)
autoplot(lm(Fv.Fm*100 ~ PDIR * Shade.Index, data = D), label.size = 3)

Forest Plot of Coefficients

• Alternative to a table for regression results!
• Can show comparisons between two competing models!
• Scale= TRUE standardizes them (so this is commonly misinterpreted to be a way of ascertaining their relative importance, but this is incorrect!)

plot_summs(Model.Final, center=FALSE, scale=TRUE, plot.distributions = TRUE)
## Loading required namespace: broom.mixed

plot_summs(Model.Final, center=TRUE, plot.distributions = TRUE, rescale.distributions = TRUE)

Compare models

Does the full model explain significantly more variation than by chance?

Plot predicted values!

#library(broom.mixed)
#effect_plot(Model.Final, interval = TRUE, plot.points = TRUE)

————————-

SUMMARY STATS TABLES

Function

• New way to enter functions after 2019 packages update:

summary(D)# means, min, max, etc for all variables
##       Time        Moss     Moss.number      Moss.sub            Microsite 
##  Min.   :0   10     : 1   Min.   : 1.00   Length:92          BB-F-37 : 1  
##  1st Qu.:0   11     : 1   1st Qu.:23.75   Class :character   BB-F-38a: 1  
##  Median :0   12     : 1   Median :48.50   Mode  :character   BB-F-39a: 1  
##  Mean   :0   13a    : 1   Mean   :48.08                      BB-F-40a: 1  
##  3rd Qu.:0   14a    : 1   3rd Qu.:72.25                      BB-F-41a: 1  
##  Max.   :0   15a    : 1   Max.   :96.00                      BB-F-42a: 1  
##              (Other):86                                      (Other) :86  
##  Old.Habitat              Habitat       Index         Elevation    Site  
##  Length:92          Interspace:17   Min.   : 3.00   Min.   : 890   1:24  
##  Class :character   Dripline  :27   1st Qu.: 8.00   1st Qu.: 890   2:34  
##  Mode  :character   Canopy    :48   Median :14.00   Median :1670   3:34  
##                                     Mean   :14.67   Mean   :1614         
##                                     3rd Qu.:20.50   3rd Qu.:2070         
##                                     Max.   :26.00   Max.   :2070         
##                                                                          
##           Veg.Type  Topography  Shade.Index    Summer.Shade.Index
##  Low-scrubland:24   S:23       Min.   : 6.00   Min.   : 2.00     
##  Mid-shrubland:34   F:35       1st Qu.:12.00   1st Qu.: 9.00     
##  High-woodland:34   N:34       Median :20.00   Median :23.00     
##                                Mean   :18.09   Mean   :21.21     
##                                3rd Qu.:23.00   3rd Qu.:32.25     
##                                Max.   :27.00   Max.   :42.00     
##                                                                  
##       PDIR            Slope            Aspect        VegPercent    
##  Min.   :0.7717   Min.   : 0.000   Min.   :  0.0   Min.   :  2.78  
##  1st Qu.:0.9036   1st Qu.: 2.000   1st Qu.: 90.0   1st Qu.: 24.30  
##  Median :0.9846   Median : 7.000   Median :200.0   Median : 38.89  
##  Mean   :0.9478   Mean   : 7.587   Mean   :193.2   Mean   : 41.46  
##  3rd Qu.:0.9997   3rd Qu.:12.000   3rd Qu.:270.0   3rd Qu.: 61.11  
##  Max.   :1.0496   Max.   :18.000   Max.   :357.0   Max.   :100.00  
##                                    NA's   :4                       
##     SolarVeg          Camera           Pair       Notes.Urgent      
##  Min.   :  4.17   Min.   : 1.00   Min.   : 1.00   Length:92         
##  1st Qu.: 29.17   1st Qu.: 6.25   1st Qu.:10.25   Class :character  
##  Median : 45.83   Median :11.50   Median :19.50   Mode  :character  
##  Mean   : 46.11   Mean   :11.73   Mean   :19.37                     
##  3rd Qu.: 63.54   3rd Qu.:16.75   3rd Qu.:28.75                     
##  Max.   :100.00   Max.   :23.00   Max.   :37.00                     
##                   NA's   :70      NA's   :54                        
##   Assay.Set             Assay            Clip             Fo        
##  Length:92          Min.   :1.000   Min.   : 1.00   Min.   : 30.00  
##  Class :character   1st Qu.:1.000   1st Qu.: 6.75   1st Qu.: 79.00  
##  Mode  :character   Median :2.000   Median :12.50   Median : 95.50  
##                     Mean   :2.467   Mean   :12.42   Mean   : 98.72  
##                     3rd Qu.:3.000   3rd Qu.:18.25   3rd Qu.:119.00  
##                     Max.   :4.000   Max.   :24.00   Max.   :174.00  
##                                                                     
##        Fm              Fv             Fv.Fm              Fs        
##  Min.   : 52.0   Min.   : 12.00   Min.   :0.1300   Min.   : 29.00  
##  1st Qu.:128.8   1st Qu.: 44.75   1st Qu.:0.3490   1st Qu.: 75.25  
##  Median :182.5   Median : 91.00   Median :0.4770   Median : 91.50  
##  Mean   :193.9   Mean   : 95.16   Mean   :0.4529   Mean   : 95.90  
##  3rd Qu.:248.2   3rd Qu.:132.00   3rd Qu.:0.5673   3rd Qu.:115.00  
##  Max.   :419.0   Max.   :265.00   Max.   :0.7370   Max.   :190.00  
##                                                                    
##       Fm.L            Fo.L             Fv.L           Fv.Fm.L      
##  Min.   : 46.0   Min.   : 17.00   Min.   : 15.00   Min.   :0.2230  
##  1st Qu.:100.8   1st Qu.: 51.75   1st Qu.: 42.00   1st Qu.:0.4047  
##  Median :136.5   Median : 70.00   Median : 61.00   Median :0.4890  
##  Mean   :141.1   Mean   : 70.58   Mean   : 70.49   Mean   :0.4877  
##  3rd Qu.:174.5   3rd Qu.: 90.00   3rd Qu.: 95.25   3rd Qu.:0.5790  
##  Max.   :319.0   Max.   :149.00   Max.   :211.00   Max.   :0.7360  
##                                                                    
##     PhiPSII             qP              qNP               NPQ          
##  Min.   :0.0430   Min.   :0.1330   Min.   :-0.8570   Min.   :-0.19100  
##  1st Qu.:0.1940   1st Qu.:0.4968   1st Qu.: 0.1580   1st Qu.: 0.08375  
##  Median :0.2685   Median :0.5820   Median : 0.4100   Median : 0.31850  
##  Mean   :0.2943   Mean   :0.5847   Mean   : 0.3634   Mean   : 0.43336  
##  3rd Qu.:0.3680   3rd Qu.:0.6783   3rd Qu.: 0.5883   3rd Qu.: 0.69700  
##  Max.   :0.6930   Max.   :0.9460   Max.   : 0.7970   Max.   : 2.54200  
##                                                                        
##    Shrubindex       Topoindex         Siteindex        Multibuffer     
##  Min.   :-1.000   Min.   :-1.0000   Min.   :-1.0000   Min.   :-2.0000  
##  1st Qu.: 0.000   1st Qu.:-0.2500   1st Qu.:-1.0000   1st Qu.: 0.0000  
##  Median : 1.000   Median : 0.0000   Median : 0.0000   Median : 1.0000  
##  Mean   : 0.337   Mean   : 0.1196   Mean   : 0.1087   Mean   : 0.5652  
##  3rd Qu.: 1.000   3rd Qu.: 1.0000   3rd Qu.: 1.0000   3rd Qu.: 1.0000  
##  Max.   : 1.000   Max.   : 1.0000   Max.   : 1.0000   Max.   : 3.0000  
##                                                                        
##    Shade.Time   
##  Min.   :21.43  
##  1st Qu.:42.86  
##  Median :71.43  
##  Mean   :64.60  
##  3rd Qu.:82.14  
##  Max.   :96.43  
## 

#Stats needed for plotting
SE<- function(x) round((sqrt(var(x)/length(x))),3)
Myfunctions = list (
  Mean = function(x) round(mean(x),3),
  Max= function(x)round(max(x),3),
  Min= function(x)round(min(x),3),
  Range= function(x) round((max(x)-min(x)),3),
  SE= function(x) round((sqrt(var(x)/length(x))),3),
  Upper= function(x) round(mean(x)+SE(x),3),
  Lower= function(x) round(mean(x)-SE(x),3),
  Var = function(x) round(var(x),3),
  SD = function(x) round(sd(x),3)
)

Tables for Sun Seeker Shade Time Boxplots

• Goal: make summary tables for SolarVeg to use in plotting and results.

#Fv/Fm Tables
#Site table
T.Site<-D %>% 
  group_by(Site) %>% 
  summarise_at(vars(Shade.Time), .funs=Myfunctions)
kable(T.Site, caption = "Shade Time Stats by Site", format="pandoc", padding=2)

Shade Time Stats by Site

Site	Mean	Max	Min	Range	SE	Upper	Lower	Var	SD
1	43.452	78.571	21.429	57.143	3.709	47.161	39.743	330.154	18.170
2	63.761	85.714	35.714	50.000	2.830	66.591	60.931	272.211	16.499
3	80.357	96.429	42.857	53.571	2.024	82.381	78.333	139.340	11.804

#Topography tables
T.Topo<-D %>% 
  group_by(Topography) %>% 
  summarise_at(vars(Shade.Time), .funs=Myfunctions)
kable(T.Topo, caption = "Shade Time Stats by Topography", format="pandoc", padding=2)

Shade Time Stats by Topography

Topography	Mean	Max	Min	Range	SE	Upper	Lower	Var	SD
S	70.652	92.857	42.857	50.000	3.238	73.890	67.414	241.137	15.529
F	58.776	89.286	21.429	67.857	3.902	62.678	54.874	532.970	23.086
N	66.492	96.429	28.571	67.857	3.631	70.123	62.861	448.373	21.175

#Microhabitat tables
T.Hab<-D %>% 
  group_by(Habitat) %>% 
  summarise_at(vars(Shade.Time), .funs=Myfunctions)
kable(T.Hab, caption = "Shade Stats by Habitat", format="pandoc", padding=2)

Shade Stats by Habitat

Habitat	Mean	Max	Min	Range	SE	Upper	Lower	Var	SD
Interspace	32.773	46.429	21.429	25.000	1.894	34.667	30.879	60.962	7.808
Dripline	59.788	89.286	28.571	60.714	3.118	62.906	56.670	262.552	16.203
Canopy	78.571	96.429	42.857	53.571	1.541	80.112	77.030	113.982	10.676

###Stress (Fv/Fm) Tables for Boxplots

#Fv/Fm Tables
#Lifezone table
T.Life.Fvfm<-D %>% group_by(Site) %>% 
  summarise_at(vars(Fv.Fm), .funs=Myfunctions)
kable(T.Life.Fvfm, caption = "Fv/Fm Stats by Lifezone Site", format="pandoc", padding=2)

Fv/Fm Stats by Lifezone Site

Site	Mean	Max	Min	Range	SE	Upper	Lower	Var	SD
1	0.433	0.571	0.183	0.388	0.021	0.454	0.412	0.011	0.104
2	0.429	0.737	0.130	0.607	0.027	0.456	0.402	0.026	0.160
3	0.491	0.686	0.144	0.542	0.028	0.519	0.463	0.026	0.163

#Topography tables
T.Topo.Fvfm<-D %>% group_by(Topography) %>% 
  summarise_at(vars(Fv.Fm), .funs=Myfunctions)
kable(T.Topo.Fvfm, caption = "Fv/Fm Stats by Topography", format="pandoc", padding=2)

Fv/Fm Stats by Topography

Topography	Mean	Max	Min	Range	SE	Upper	Lower	Var	SD
S	0.432	0.656	0.20	0.456	0.026	0.458	0.406	0.016	0.125
F	0.383	0.647	0.13	0.517	0.027	0.410	0.356	0.026	0.162
N	0.539	0.737	0.25	0.487	0.018	0.557	0.521	0.012	0.107

#Microhabitat tables
T.Hab.Fvfm<-D %>% group_by(Habitat) %>% 
  summarise_at(vars(Fv.Fm), .funs=Myfunctions)
kable(T.Hab.Fvfm, caption = "Fv/Fm Stats by Habitat", format="pandoc", padding=2)

Fv/Fm Stats by Habitat

Habitat	Mean	Max	Min	Range	SE	Upper	Lower	Var	SD
Interspace	0.335	0.508	0.160	0.348	0.025	0.360	0.310	0.011	0.103
Dripline	0.414	0.647	0.130	0.517	0.028	0.442	0.386	0.021	0.144
Canopy	0.516	0.737	0.174	0.563	0.020	0.536	0.496	0.018	0.136

#—————————————– # FINAL ANOVA PLOTS ~ Fv/Fm * n= 92 (stats will change) * plot Welch’s ANVOA results on graphs? (But must choose 1 of 3 conservative parametric pairwise comparison tests.) *

Plot A - Life zone

• Graphics note: points are only jittered on the x-axis using position_dodge2().

#Veg ~ life zone boxplots

PlotA<-ggplot(D, aes(x=(Site), y=Shade.Time))+
  geom_boxplot(outlier.shape=NA, aes(x=factor(Site), fill= (Site)))+
  geom_point(aes(fill=(Site)), position=position_dodge2(width = 0.4), size=.5, 
             color="black")+
  stat_summary(fun=mean, geom="point", shape=20, size=3, color="white")+
  #geom_point(data = T.Site, aes(x=Site, y=Mean), color="white")+
  geom_errorbar(data=T.Site, aes(x=Site, y= Mean, ymin=Lower, ymax=Upper), width=.2, size=0.7, color="white") +
  ggtitle("Macroscale habitat")+
  labs(y='Shade time (% of year)', x="Elevation (life zone site)")+
  #scale_x_discrete(labels=c("1" = "Low-scrubland", "2" = "Mid-shrubland", "3"= "High-woodland"))+
  scale_x_discrete(labels=c("1" = "Low", "2" = "Mid", "3"= "High"))+
  scale_y_continuous(name = "Annual shade time (% of year)",breaks = seq(0, 100, 20), limits=c(0, 100.1), expand=c(0,0))+
   scale_fill_manual(values =c("#E6AB02", "firebrick3", "dodgerblue3"))+
  theme_linedraw()+#must be BEFORE the themes below!
  theme(legend.position = "none", text=element_text(size=18), 
        axis.title.x = element_blank()) # panel.grid.minor = element_blank()
PlotA

Plot B - Topo

 PlotB<-ggplot(D)+
    aes(x=Topography, y = Shade.Time, fill = Topography) +
    geom_boxplot(outlier.shape=NA)+ #hides outliers so they are not double-plotted
    geom_point(position = position_dodge2( width=0.4), size=.5, color="black")+
    stat_summary(fun=mean, geom="point", shape=20, size=3, color="white")+
    geom_errorbar(data=T.Topo, aes(x=Topography, y= Mean, ymin=Lower, ymax=Upper), 
                width=.2, size=0.7, color="white") +
    scale_y_continuous(name = " ",breaks = seq(0, 100, 20), 
                        limits=c(0, 100.1), expand=c(0,0))+
    scale_x_discrete(name = "Topography zone", 
                     labels=c("S-facing", "Flat", "N-facing")) +
   ggtitle("Mesoscale habitat") +
      theme_linedraw() +
      theme(text = element_text(size=18), legend.position = "none") +
         scale_fill_manual(values = c("#fe9929", "#666666","darkorchid"))+
       theme(panel.background = element_rect(fill = "white"), 
             axis.title.x = element_blank())# panel.grid.minor = 
             # element_blank())#axis.title.x = element_blank(),
PlotB

Plot C - Microhabitat

Microhabitat shade time

PlotC<-ggplot(D)+
    aes(x=Habitat, y =Shade.Time, fill = Habitat) +
    geom_boxplot(outlier.shape=NA)+
    geom_point(position = position_dodge2(width=0.4), size=0.5, color="black")+
    stat_summary(fun=mean, geom="point", shape=20, size=3, color="white")+
    geom_errorbar(data=T.Hab, aes(x=Habitat, y= Mean, ymin=Lower, ymax=Upper),
                width=.15, size=0.7, color= "white") +
    scale_y_continuous(name = " ",breaks = seq(0, 100, 20), 
                       limits=c(0, 100.1), expand=c(0,0))+
    scale_x_discrete(name = "Microhabitat type") +
    ggtitle("Microscale habitat", )+
    theme_linedraw() +
    theme(text = element_text(size=18), legend.position = "none") +
       scale_fill_manual(values = c("tan", "yellowgreen","darkgreen"))+
    theme(panel.background = element_rect(fill = "white"), 
          axis.title.x = element_blank())#panel.grid.minor = 
          #element_blank())#
PlotC

Plot D - Stress vs Lifezone

• Outlier (Phi-PSII) is Moss 34 at Site 2. How could this reading be SO high? Was Fv/Fm high also? Yes: 0.737.
• Question: what dimensions to save for publication? (fig.width=4, fig.height = 1.7)

#fig.width = 3.5, fig.height = 3,
#Fv/Fm life zone boxplots
PlotD<-ggplot(D, aes(x=Site, y=Fv.Fm, fill = Site))+
  geom_boxplot(outlier.shape=NA)+ 
  geom_point(position = position_dodge2(width=0.4),size=0.5,color="black")+
  #stat_summary(fun.y=mean, geom="point", shape=20, size=3, color="white")+
  geom_point(data = T.Life.Fvfm, aes(x=Site, y=Mean), color="white")+
  geom_errorbar(data=T.Life.Fvfm, aes(x=Site, y= Mean, ymin=Lower, ymax=Upper), 
                width=.2, size=0.7, color="white") +
  labs(y='\u2190 Stress (Fv/Fm)', x="Elevation-life zone (site)")+
  scale_x_discrete(labels=c("1" = "Low", "2" = "Mid", "3"= "High"))+
  scale_y_continuous(breaks = seq(0, 0.78, 0.1),expand=c(0,0), limits=c(0,0.76))+
  scale_fill_manual(values =c("#E6AB02", "firebrick3", "dodgerblue3"), guide="none")+
  theme_linedraw()+# panel.grid.major = element_blank(),
  theme(legend.position = "none", text=element_text(size=18), panel.grid.minor = 
         element_blank())#axis.title.x = element_blank()
  
PlotD

Plot E - Stress vs. Topo

#Fv/Fm topo boxplots
PlotE<-ggplot(D)+
    aes(x=Topography, y = Fv.Fm, fill = Topography) +
    geom_boxplot(outlier.shape=NA)+
    geom_point(position = position_dodge2( width=0.4), size=0.5, color="black")+
    stat_summary(fun=mean, geom="point", shape=20, size=3, color="white")+
    geom_errorbar(data=T.Topo.Fvfm, aes(x=Topography, y= Mean, ymin=Lower, ymax=Upper), width=.2, size=0.7, color="white") +
    scale_y_continuous(name = "",
       breaks = seq(0, 0.7, 0.1), expand=c(0,0), limits=c(0,0.76)) +
    scale_x_discrete(name = "Topography zone", 
                     labels=c("S-facing", "Flat", "N-facing")) +
   #ggtitle("Summer stress assay (Series 1.1)") +
      theme_linedraw() +
      theme(text = element_text(size=18), legend.position = "none") +
         scale_fill_manual(values = c("#fe9929", "#666666","darkorchid"))+
       theme(panel.background = element_rect(fill = "white"), panel.grid.minor = 
              element_blank())#axis.title.x = element_blank(),
PlotE# Y axis: "\u2190 Stress (Fv/Fm)"

Plot F - Stress vs. Microhabitat

# Welch One way ANOVA test
AovF <- D %>% welch_anova_test(Fv.Fm ~ Habitat)
AovF
## # A tibble: 1 × 7
##   .y.       n statistic   DFn   DFd          p method     
## * <chr> <int>     <dbl> <dbl> <dbl>      <dbl> <chr>      
## 1 Fv.Fm    92      16.6     2  44.5 0.00000414 Welch ANOVA
# Pairwise comparisons (Games-Howell)
 
PWC<- D %>% games_howell_test(Fv.Fm ~ Habitat, detailed=TRUE)

PlotF<-ggplot(D)+
    aes(x=Habitat, y =Fv.Fm, fill = Habitat) +
    geom_boxplot(outlier.shape=NA)+
    geom_point(position = position_dodge2( width=0.4), size=.5, color="black")+
    stat_summary(fun=mean, geom="point", shape=20, size=3, color="white")+
    geom_errorbar(data=T.Hab.Fvfm, aes(x=Habitat, y= Mean, ymin=Lower, ymax=Upper), 
                width=.15, size=0.7, color= "white") +
    scale_y_continuous(name = " ",breaks = seq(0, 0.75, 0.1),
                       expand=c(0,0), limits=c(0,0.76)) +
    scale_x_discrete(name = "Microhabitat type") +
  stat_pwc(method = "games_howell_test")+
    labs(subtitle = get_test_label(AovF, detailed = TRUE),
    caption = get_pwc_label(PWC))+
  # ggtitle("Micro-habitat buffering") +
      theme_linedraw() +
      theme(text = element_text(size=18), legend.position = "none") +
         scale_fill_manual(values = c("tan", "yellowgreen","darkgreen"))+
       theme(panel.background = element_rect(fill = "white"), 
             panel.grid.minor = element_blank())#axis.title.x = element_blank()
##Stress text with arrow: \u2190 Stress (Fv/Fm) 
#PlotF + stat_pvalue_manual(PWC, hide.ns = TRUE) 

FINAL FIGURE ANOVA boxplots: Shade & Stress vs 3 scales of habitat structure

• updated with n= 92

#Plot A-F together
PlotsA_F <-ggarrange(PlotA, PlotB, PlotC, PlotD, PlotE, PlotF, labels=c("A","B","C","D","E", "F"), widths=1:1, ncol=3, nrow=2, font.label= list(size=18)) #Add letters & up the size
#PlotsA_F<- annotate_figure(PlotsA_F, top = text_grob("Shade time & summer stress vs. 3 scales of habitat buffering", size = 16, rot=0, just= "center")) #Add bottom/top titles
PlotsA_F

#——————————————- # VARIANCE TESTS

• Goal: Test null of equal variance across factor levels for shade & stress.

• Test: Fligner-Killeen homogeneity of variance test: stats::fligner.test(). For when groups lack normal distributions.

• Results:

  • Only Fv/Fm by Topography is heteroskedastic
  • Fv/Fm variance was equal for Life zone & Habitat.
  • Shade variances are all heteroskedastic

#FvFm variance tests for life zone sites (Veg.Type)

fligner.test(Fv.Fm ~ Veg.Type, data = D) #Fligner-Killeen:med chi-squared = chi-squared = 4.7903, df = 2, p-value = 0.09116 #Equal!

fligner.test(Fv.Fm ~ Topography, data = D) #median chi-squared = 7.5076, df = 2, p-value = 0.02343 #Not equal

fligner.test(Fv.Fm ~ Habitat, data = D)#med chi-squared = 1.3133, df = 2, p-value = 0.5186 #equal

#Shade.time tests (Shade.Index/100)
fligner.test(Shade.Time ~ Veg.Type, data = D) #Fligner-Killeen:med chi-squared = chi-squared = 8.2128, df = 2, p-value = 0.01647 #Unequal

fligner.test(Shade.Time ~ Topography, data = D) #median chi-squared  6.4838, df = 2, p-value = 0.03909 #Unequal

fligner.test(Shade.Time ~ Habitat, data = D)#med chi-squared == 8.6992, df = 2, p-value = 0.01291 #unequal
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  Fv.Fm by Veg.Type
## Fligner-Killeen:med chi-squared = 4.7903, df = 2, p-value = 0.09116
## 
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  Fv.Fm by Topography
## Fligner-Killeen:med chi-squared = 7.5076, df = 2, p-value = 0.02343
## 
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  Fv.Fm by Habitat
## Fligner-Killeen:med chi-squared = 1.3133, df = 2, p-value = 0.5186
## 
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  Shade.Time by Veg.Type
## Fligner-Killeen:med chi-squared = 8.2128, df = 2, p-value = 0.01647
## 
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  Shade.Time by Topography
## Fligner-Killeen:med chi-squared = 6.4838, df = 2, p-value = 0.03909
## 
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  Shade.Time by Habitat
## Fligner-Killeen:med chi-squared = 8.6992, df = 2, p-value = 0.01291

#——————————————- # WELCH’s ANOVAs

• Report all Welch’s ANOVAs which perform best for unequal sample sizes by group level and unequal variance (if present or not).

• use stats::oneway.test() with argument var.equal = FALSE, for Welch Correction, an approximate method of Welch (1951), which generalizes the 2-sample Welch test to the case of arbitrarily many samples. Welch’s correction reduces the denominator degrees of freedom in the F-test as a penalty for unequal variance. I can’t find published documentation on the mechanics of this test or how the Degrees of Freedom are calculated!

Fv.Fm vs Life zone

Welch’s Anova - use this Life zone model, which adjusts for unequal variance even though variance was equal, this test is best for unequal sample sizes by life zone.

library(purrr)
#Welch's ANOVA
W.Aovs<- D %>% select(Fv.Fm) %>% #select variables
  purrr::map_df(~broom::tidy(oneway.test(. ~ Veg.Type, var.equal= FALSE, data= D)),.id="Variable")#run classic anova
## Multiple parameters; naming those columns num.df, den.df
W.Aovs<-W.Aovs %>% select(-method)

#Extract Cohen's F from batch
CF<-D %>% select(Fv.Fm) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Veg.Type, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$cohens.f
  })
CF<-as.data.frame(CF) %>% round(., 3)#make df
CF<-t(CF)#transpose
CF<- CF[,1]
W.Aovs<-W.Aovs %>% mutate(Cohens_F = CF)#Add to output

#Extract eta squared (R-squared)
R.sq<-D %>% select(Fv.Fm) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Veg.Type, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$etasq
  })
R.sq<-as.data.frame(R.sq) %>% round(., 4)#make dataframe
R.sq<-t(R.sq)#transpose
W.Aovs<-W.Aovs %>% mutate(Rsq = R.sq)#Add to output

#Print Pretty Table
kable(W.Aovs, caption = "Welch's ANOVAs Fv/Fm vs Life zone", digits = c(1, 1, 3, 3, 4, 5, 3)) %>% kable_styling()# table.attr = "style = \"color: black;\""

Welch’s ANOVAs Fv/Fm vs Life zone

Variable	num.df	den.df	statistic	p.value	Cohens_F	Rsq
Fv.Fm	2	59.284	1.633	0.2039	0.198	0.038

Fv.Fm ~ Topo

#Welch's ANOVA - use var.equal = FALSE
W.Aovs<- D %>% select(Fv.Fm) %>% #select variables
  purrr::map_df(~broom::tidy(oneway.test(. ~ Topography, var.equal= FALSE, data= D)),.id="Variable")#run classic anova
## Multiple parameters; naming those columns num.df, den.df
W.Aovs<-W.Aovs %>% select(-method)

#Extract Cohen's F from batch
CF<-D %>% select(Fv.Fm) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Topography, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$cohens.f
  })
CF<-as.data.frame(CF) %>% round(., 3)#make df
CF<-t(CF)#transpose
CF<- CF[,1]
W.Aovs<-W.Aovs %>% mutate(Cohens_F = CF)#Add to output

#Extract eta squared (R-squared)
R.sq<-D %>% select(Fv.Fm) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Topography, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$etasq
  })
R.sq<-as.data.frame(R.sq) %>% round(., 4)#make dataframe
R.sq<-t(R.sq)#transpose
W.Aovs<-W.Aovs %>% mutate(Rsq = R.sq)#Add to output

#Print Pretty Table
kable(W.Aovs, caption = "Welch's ANOVAs Topography vs Fv/Fm", digits = c(1, 1, 3, 3, 4, 5, 3)) %>% kable_styling()# table.attr = "style = \"color: black;\""

Welch’s ANOVAs Topography vs Fv/Fm

Variable	num.df	den.df	statistic	p.value	Cohens_F	Rsq
Fv.Fm	2	53.682	12.759	0	0.515	0.21

Fv.Fm ~ Habitat

#1-way Welch's ANOVA - use var.equal = FALSE
W.Aovs<- D %>% select(Fv.Fm) %>% #select variables
  purrr::map_df(~broom::tidy(oneway.test(. ~ Habitat, var.equal= FALSE, data= D)),.id="Variable")#run classic anova
## Multiple parameters; naming those columns num.df, den.df
W.Aovs<-W.Aovs %>% select(-method)

#Extract Cohen's F from batch
CF<-D %>% select(Fv.Fm) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Habitat, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$cohens.f
  })
CF<-as.data.frame(CF) %>% round(., 3)#make df
CF<-t(CF)#transpose
CF<- CF[,1]
W.Aovs<-W.Aovs %>% mutate(Cohens_F = CF)#Add to output

#Extract eta squared (R-squared)
R.sq<-D %>% select(Fv.Fm) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Habitat, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$etasq
  })
R.sq<-as.data.frame(R.sq) %>% round(., 4)#make dataframe
R.sq<-t(R.sq)#transpose
W.Aovs<-W.Aovs %>% mutate(Rsq = R.sq)#Add to output

#Print Pretty Table
kable(W.Aovs, caption = "Welch's ANOVAs Habitat vs Fv/Fm", digits = c(1, 1, 3, 3, 4, 5, 3)) %>% kable_styling()#table.attr = "style = \"color: black;\""

Welch’s ANOVAs Habitat vs Fv/Fm

Variable	num.df	den.df	statistic	p.value	Cohens_F	Rsq
Fv.Fm	2	44.537	16.584	0	0.547	0.23

Games-Howell multiple comparisons

Results * S vs F Topography zones not sign. different * Interspace vs Dripline n.s.

#Multiple comparisons - Games-Howell 
#library(rstatix)
games_howell_test(data=D, formula= Fv.Fm ~ Topography, conf.level=0.95, detailed=F)
## # A tibble: 3 × 8
##   .y.   group1 group2 estimate conf.low conf.high     p.adj p.adj.signif
## * <chr> <chr>  <chr>     <dbl>    <dbl>     <dbl>     <dbl> <chr>       
## 1 Fv.Fm S      F       -0.0484  -0.139     0.0425 0.411     ns          
## 2 Fv.Fm S      N        0.107    0.0294    0.184  0.005     **          
## 3 Fv.Fm F      N        0.155    0.0758    0.235  0.0000471 ****

games_howell_test(data=D, formula= Fv.Fm ~ Habitat, conf.level=0.95, detailed=F)
## # A tibble: 3 × 8
##   .y.   group1     group2   estimate conf.low conf.high      p.adj p.adj.signif
## * <chr> <chr>      <chr>       <dbl>    <dbl>     <dbl>      <dbl> <chr>       
## 1 Fv.Fm Interspace Dripline   0.0789  -0.0120     0.170 0.1        ns          
## 2 Fv.Fm Interspace Canopy     0.181    0.104      0.259 0.00000476 ****        
## 3 Fv.Fm Dripline   Canopy     0.103    0.0207     0.184 0.011      *

Shade vs Life zone

Welch’s Anova - use this Life zone model, which adjusts for unequal variance because variance was unequal & this test is best for unequal sample sizes by life zone.

#library(purrr)
#Welch's ANOVAs - use var.equal = FALSE
W.Aovs<- D %>% select(Shade.Time) %>% #select variables
  purrr::map_df(~broom::tidy(oneway.test(. ~ Veg.Type, var.equal= FALSE, data= D)),.id="Variable")#run classic anova
## Multiple parameters; naming those columns num.df, den.df
W.Aovs<-W.Aovs %>% select(-method)

#Extract Cohen's F from batch
CF<-D %>% select(Shade.Time) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Veg.Type, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$cohens.f
  })
CF<-as.data.frame(CF) %>% round(., 3)#make df
CF<-t(CF)#transpose
CF<- CF[,1]
W.Aovs<-W.Aovs %>% mutate(Cohens_F = CF)#Add to output

#Extract eta squared (R-squared)
R.sq<-D %>% select(Shade.Time) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Veg.Type, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$etasq
  })
R.sq<-as.data.frame(R.sq) %>% round(., 4)#make dataframe
R.sq<-t(R.sq)#transpose
W.Aovs<-W.Aovs %>% mutate(Rsq = R.sq)#Add to output

#Print Pretty Table
kable(W.Aovs, caption = "Welch's ANOVAs Shade vs Life zone", digits = c(1, 1, 3, 3, 4, 5, 3)) %>% kable_styling()#table.attr = "style = \"color: black;\""

Welch’s ANOVAs Shade vs Life zone

Variable	num.df	den.df	statistic	p.value	Cohens_F	Rsq
Shade.Time	2	51.287	40.57	0	0.952	0.476

Shade vs Topo

library(purrr)
#Welch's ANOVA - use var.equal = FALSE
W.Aovs<- D %>% select(Shade.Time) %>% #select variables
  purrr::map_df(~broom::tidy(oneway.test(. ~ Topography, var.equal= FALSE, data= D)),.id="Variable")#run classic anova
## Multiple parameters; naming those columns num.df, den.df
W.Aovs<-W.Aovs %>% select(-method)

#Extract Cohen's F from batch
CF<-D %>% select(Shade.Time) %>% 
  map_dfr(., function(x) {
     Mod1<- aov(x ~ Topography, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$cohens.f
  })
CF<-as.data.frame(CF) %>% round(., 3)#make df
CF<-t(CF)#transpose
CF<- CF[,1]
W.Aovs<-W.Aovs %>% mutate(Cohens_F = CF)#Add to output

#Extract eta squared (R-squared)
R.sq<-D %>% select(Shade.Time) %>% 
  map_dfr(., function(x) {
     Mod1<- aov(x ~ Topography, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$etasq
  })
R.sq<-as.data.frame(R.sq) %>% round(., 4)#make dataframe
R.sq<-t(R.sq)#transpose
W.Aovs<-W.Aovs %>% mutate(Rsq = R.sq)#Add to output

#Print Pretty Table
kable(W.Aovs, caption = "Welch's ANOVAs Topography vs Shade", digits = c(1, 1, 3, 3, 4, 5, 3)) %>% kable_styling()#table.attr = "style = \"color: black;\""

Welch’s ANOVAs Topography vs Shade

Variable	num.df	den.df	statistic	p.value	Cohens_F	Rsq
Shade.Time	2	58.491	2.724	0.0739	0.237	0.053

Shade ~ Habitat

#1-way Welch's ANOVA - use var.equal = FALSE
W.Aovs<- D %>% select(Shade.Time) %>% #select variables
  purrr::map_df(~broom::tidy(oneway.test(. ~ Habitat, var.equal= FALSE, data= D)),.id="Variable")#run classic anova
## Multiple parameters; naming those columns num.df, den.df
W.Aovs<-W.Aovs %>% select(-method)

#Extract Cohen's F from batch
CF<-D %>% select(Shade.Time) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Habitat, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$cohens.f
  })
CF<-as.data.frame(CF) %>% round(., 3)#make df
CF<-t(CF)#transpose
CF<- CF[,1]
W.Aovs<-W.Aovs %>% mutate(Cohens_F = CF)#Add to output

#Extract eta squared (R-squared)
R.sq<-D %>% select(Shade.Time) %>% 
  map_df(., function(x) {
     Mod1<- aov(x ~ Habitat, data=D)
     Out<-sjstats::anova_stats(Mod1)#car pagackage - calculates Cohen's F
     Out$etasq
  })
R.sq<-as.data.frame(R.sq) %>% round(., 4)#make dataframe
R.sq<-t(R.sq)#transpose
W.Aovs<-W.Aovs %>% mutate(Rsq = R.sq)#Add to output

#Print Pretty Table
kable(W.Aovs, caption = "Welch's ANOVAs Habitat vs Shade", digits = c(1, 1, 3, 3, 4, 5, 3)) %>% kable_styling()#table.attr = "style = \"color: black;\""

Welch’s ANOVAs Habitat vs Shade

Variable	num.df	den.df	statistic	p.value	Cohens_F	Rsq
Shade.Time	2	43.953	173.329	0	1.438	0.674

Games-Howell multiple comparisons

Results * *

#Multiple comparisons - Games-Howell 
#library(rstatix)
games_howell_test(data=D, formula= Shade.Time ~ Veg.Type, conf.level=0.95, detailed=F)
## # A tibble: 3 × 8
##   .y.        group1     group2 estimate conf.low conf.high    p.adj p.adj.signif
## * <chr>      <chr>      <chr>     <dbl>    <dbl>     <dbl>    <dbl> <chr>       
## 1 Shade.Time Low-scrub… Mid-s…     20.3     9.01      31.6 2.12e- 4 ***         
## 2 Shade.Time Low-scrub… High-…     36.9    26.6       47.2 5.39e-10 ****        
## 3 Shade.Time Mid-shrub… High-…     16.6     8.23      25.0 3.61e- 5 ****

games_howell_test(data=D, formula= Shade.Time ~ Habitat, conf.level=0.95, detailed=F)
## # A tibble: 3 × 8
##   .y.        group1     group2  estimate conf.low conf.high   p.adj p.adj.signif
## * <chr>      <chr>      <chr>      <dbl>    <dbl>     <dbl>   <dbl> <chr>       
## 1 Shade.Time Interspace Dripli…     27.0     18.1      35.9 1.56e-8 ****        
## 2 Shade.Time Interspace Canopy      45.8     39.8      51.7 0       ****        
## 3 Shade.Time Dripline   Canopy      18.8     10.3      27.3 1.04e-5 ****

#————————– # ANCOVA * Updated with n=92

Interaction Plot

#Fv/Fm plot
ggplot(D, aes(Shade.Time, Fv.Fm, color= (Topography)))+
    facet_grid(factor(Veg.Type, levels=c("High-woodland", "Mid-shrubland", "Low-scrubland"))~.)+
     geom_smooth(method=lm, fill= "grey", formula = y ~ x)+
    geom_point(aes(color=Topography), size= 2) +
    labs(x = "Annual shade time (% of the year)", 
         y="\u2190 Moss stress (Fv/Fm)")+
    scale_color_manual(values = c("#fe9929", "#666666","darkorchid"),
        name="Topography", labels=c("S-facing","Flat", "N-facing"))+
    scale_y_continuous(breaks = seq(0,0.9,by = 0.2), limits=c(0,0.9), 
                       expand=c(0,0),minor_breaks = NULL)+
    scale_x_continuous(breaks=seq(0,100, by= 10), minor_breaks = waiver())+
    theme_linedraw()+
   #legend.text = element_text(size=15), text=element_text(size = 10), 
    theme(legend.position = "top",  legend.background = element_blank(),  
    panel.background = element_rect(fill = "antiquewhite2"),
             panel.grid.minor =  element_blank(), text=element_text(size = 15))

# panel.grid.major = element_blank(),

Center covariate variable

• Just as in regression, we need to center the continuous predictor variable (Shade.Time) to remove structural multicollinearity resulting from a significant interaction term!
• Can center with the scale(X, center = TRUE, scale = FALSE), but I’m not sure how to combine with mutate()

D<-D %>% mutate(Shade.Time.C= Shade.Time - mean(Shade.Time, na.rm = T))

SITE 1 MODEL Fv/Fm

• Note - table output is done with pander() at present. Replace with kable after updating MAcOS.
• 2023 QUESTION: if we are doing ANCOVA, why is step 1 to fit a linear model (lm)?? I guess because the actual ANCOVA functions (aov(), Anova(), oneway.test) can take input of the lm, so rather than type it a bunch of times, assign it to an object once.
• But I’m not sure what Sums of Squares (SS) Anova to use to determine significance of the factors - type I, II, III?
• I’m also confused how to perform post-hoc tests with a “HCCM-corrected matrix”
• I deleted the Welch’s ANOVA which is ONLY for 1-way Anova, no covariates. The only output from Welch’s ANOVA is an adjusted F-stat & DF, so how to get the table?
• the final reduced model is called Mod.Topo (Shade.Time covariate was not significant, so it’s reduced to an intercepts-only regression (anova) model.

Type I Errors ANCOVA & OLS Regression

• Use only Site 1 data (Low-scrubland) for this model = D1
• Variance tests say Topozones are equal (ANCOVA assumption).
• Use centered Shade.Time.C so you avoid structural multicollinearity.

#Site 1 - Variance Test for ANOVA (grouped factor data)
D1<-subset(D, Site ==1)#Subset Site 1 
boxplot((Fv.Fm) ~ Topography, data = D1)

fligner.test(Fv.Fm ~ Topography, data= D1)# equal
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  Fv.Fm by Topography
## Fligner-Killeen:med chi-squared = 1.7132, df = 1, p-value = 0.1906
leveneTest(data = D1, log(Fv.Fm) ~ Topography)#equal medians test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  1  3.7964 0.06423 .
##       22                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#BP Breuch-Pagan Test for heteroskedasticity (in regression)
#bptest(Mod1)#not heteroskedastic but close P = 0.091

#Type I ANCOVA - 
Mod1 <- lm(Fv.Fm ~ Topography*Shade.Time.C, data = D1)
#kable(Mod1, caption = "Site 1: Type I ANCOVA Moss Stress vs Topo & Shade time", digits=c(1, 1, 4, 4, 2, 4))
summary.lm(Mod1)#Get R2
## 
## Call:
## lm(formula = Fv.Fm ~ Topography * Shade.Time.C, data = D1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.167821 -0.037291 -0.006309  0.042979  0.135179 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               0.388575   0.059630   6.516 2.37e-06 ***
## TopographyN               0.114869   0.064573   1.779   0.0905 .  
## Shade.Time.C              0.001048   0.001770   0.592   0.5605    
## TopographyN:Shade.Time.C -0.001636   0.002168  -0.755   0.4593    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07216 on 20 degrees of freedom
## Multiple R-squared:  0.5832, Adjusted R-squared:  0.5207 
## F-statistic: 9.329 on 3 and 20 DF,  p-value: 0.0004609
summary.aov(Mod1)#Type I ANOVA table
##                         Df  Sum Sq Mean Sq F value Pr(>F)    
## Topography               1 0.14276 0.14276  27.417  4e-05 ***
## Shade.Time.C             1 0.00001 0.00001   0.002  0.967    
## Topography:Shade.Time.C  1 0.00297 0.00297   0.569  0.459    
## Residuals               20 0.10414 0.00521                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R2_site1<-0.5832 #(shade.index was 0.3908)
plot(Mod1)# Looks good except Moss 3 has (high?) leverage = 0.5 

### Type I ANOVA Full Model

anova(Mod1)#Type I ANOVA table
## Analysis of Variance Table
## 
## Response: Fv.Fm
##                         Df   Sum Sq  Mean Sq F value    Pr(>F)    
## Topography               1 0.142758 0.142758 27.4170 4.004e-05 ***
## Shade.Time.C             1 0.000009 0.000009  0.0018    0.9667    
## Topography:Shade.Time.C  1 0.002965 0.002965  0.5695    0.4593    
## Residuals               20 0.104139 0.005207                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Type II Robust ANOVA

• Step 2 - Compare to Robust Type II ANOVA
• use if interaction is not important (or almost always according to car::Anova help file!)
• No intercept in Type II.
• Same F-stat & P-value as Type III.
• Which contrasts to do? (Topography levels don’t need intercept model)
• Just report regression line equations for each line from summary()

Anova(Mod1, type="II", white.adjust = "hc3")
## Coefficient covariances computed by hccm()
## Analysis of Deviance Table (Type II tests)
## 
## Response: Fv.Fm
##                         Df       F    Pr(>F)    
## Topography               1 23.0000 0.0001101 ***
## Shade.Time.C             1  0.1683 0.6859503    
## Topography:Shade.Time.C  1  0.3191 0.5784437    
## Residuals               20                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#tidy(AII<-Anova(Mod1, white.adjust=TRUE, type="2")), #KEY FIX!!
     # caption = "Site 1: White-Huber (HCCM) Type II SS T0-stress vs Topo & Shade" 
     #must adjust contrasts or result is nonsensical with intercept term!

Regression of significant factors only

• Topography has equal var: don’t need to use White.adjust HCCM-ANCOVA but Welch’s ANOVA is best for unbalanced designs (Topo has different sample sizes/level.)
  
• HCCM = (heteroskedasticity-corrected covariance matrix)
  
• adjusts the F-stat & P-value, no other SS stuff changes.
• FIX: kable tables are NOT working with this older Mac version. Update and try.
• I’m using pander() table output for now. Titles appear at the bottom of tables.

#Step 1 - Build an aov or lm model first (not an anova() object!)
Mod1b<-lm(data= D1, Fv.Fm ~ Topography)#Reduced ANOVA model
summary.lm(Mod1b)#regression model intercept is the dummy variable: Topo= Flat zone is the intercept = 0.3555= mean Fv/Fm for this factor level.
## 
## Call:
## lm(formula = Fv.Fm ~ Topography, data = D1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.172500 -0.030813  0.000375  0.034688  0.130500 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.35550    0.02014  17.649 1.79e-14 ***
## TopographyN  0.15425    0.02849   5.415 1.94e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06978 on 22 degrees of freedom
## Multiple R-squared:  0.5713, Adjusted R-squared:  0.5518 
## F-statistic: 29.32 on 1 and 22 DF,  p-value: 1.941e-05
plot(Mod1b)

Step 3 - Robust Type II ANOVA

• use if interaction is not important (or almost always according to car::Anova help file!)
• No intercept in Type II.
• Same F-stat & P-value as Type III.
• Which contrasts to do? (Topography levels don’t need intercept model)
• Just report regression line equations for each line from summary()

Mod1a<-lm(data= D1, Fv.Fm ~ Topography+Shade.Time.C)
Anova(Mod1a, type="II", white.adjust = "hc3")
## Coefficient covariances computed by hccm()
## Analysis of Deviance Table (Type II tests)
## 
## Response: Fv.Fm
##              Df       F    Pr(>F)    
## Topography    1 21.5680 0.0001394 ***
## Shade.Time.C  1  0.0019 0.9654415    
## Residuals    21                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#tidy(AII<-Anova(Mod1, white.adjust=TRUE, type="2")), #KEY FIX!!
     # caption = "Site 1: White-Huber (HCCM) Type II SS T0-stress vs Topo & Shade" 
     #must adjust contrasts or result is nonsensical with intercept term!

ROBUST HCCM with ANOVA (Type II)

• Redo with Anova(Mod, white.adjust = “hc3”, type = II)
• HCCM changes Anova() output to a “Analysis of Deviance” table with no SS’s…?
• Not possible to get the robust ANOVA table…?
• Have to save the Varcov matrix as an object using car::hccm()
• Have to use different function to get the coefficients: coeftest(Mod, varc).

### Robust additive Model with shade
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
#save VarCov matric
Mod1<-lm(data= D1, Fv.Fm ~ Topography+Shade.Time.C)
Cov1 <- hccm(Mod1, type="hc3") #needs package 'car'
Mod1.HC3 <- coeftest(Mod1, vcov.=Cov1)#needs lmtest package 
tidy(Mod1.HC3)#t tests of HCCM model coefficients
## # A tibble: 3 × 5
##   term           estimate std.error statistic      p.value
##   <chr>             <dbl>     <dbl>     <dbl>        <dbl>
## 1 (Intercept)   0.354      0.0391      9.05   0.0000000108
## 2 TopographyN   0.155      0.0334      4.64   0.000139    
## 3 Shade.Time.C -0.0000431  0.000984   -0.0438 0.965

#Run Robust ANCOVA
Anova(Mod1, white.adjust="hc3", type="II")#HCCM type HC3 Analysis of Deviance Table
## Coefficient covariances computed by hccm()
## Analysis of Deviance Table (Type II tests)
## 
## Response: Fv.Fm
##              Df       F    Pr(>F)    
## Topography    1 21.5680 0.0001394 ***
## Shade.Time.C  1  0.0019 0.9654415    
## Residuals    21                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Mod.Topo<-Anova(Mod1, white.adjust="hc3", type="II")#HCCM type HC3
## Coefficient covariances computed by hccm()
Mod.Topo# No SumSqs table only F-tests
## Analysis of Deviance Table (Type II tests)
## 
## Response: Fv.Fm
##              Df       F    Pr(>F)    
## Topography    1 21.5680 0.0001394 ***
## Shade.Time.C  1  0.0019 0.9654415    
## Residuals    21                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#summary.lm(Mod.Topo)#Rsq= 0.5713 plots the regression output
#summary.aov(Mod.Topo)#ANOVA Table (but does white.adjust change these SS's??)

pander(Anova(Mod1, white.adjust="hc3", type="II"), caption = "Site 1: Type II HCCM: Topography x shade")
## Coefficient covariances computed by hccm()

Site 1: Type II HCCM: Topography x shade ### Reduced Model R-squared * Fit reduced model & test if interaction improves fit significantly - No!

	Df	F	Pr(>F)
Topography	1	21.57	0.0001394
Shade.Time.C	1	0.001922	0.9654
Residuals	21	NA	NA

Mod1.H0<-lm(Fv.Fm~ 1, data = D1)
Mod.1red<-lm(Fv.Fm ~Topography, data = D1)
Cov1 <- hccm(Mod.1red, type="hc3") #needs package 'car'
summary.lm(Mod.1red)
## 
## Call:
## lm(formula = Fv.Fm ~ Topography, data = D1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.172500 -0.030813  0.000375  0.034688  0.130500 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.35550    0.02014  17.649 1.79e-14 ***
## TopographyN  0.15425    0.02849   5.415 1.94e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06978 on 22 degrees of freedom
## Multiple R-squared:  0.5713, Adjusted R-squared:  0.5518 
## F-statistic: 29.32 on 1 and 22 DF,  p-value: 1.941e-05
#anova(Mod.1red, Mod.1)#Proof of concept, does not improve fit
Cov1 <- hccm(Mod.1red, type="hc3") #needs package 'car'
waldtest(Mod1.H0, Mod.1red, vcov=Cov1)#Overall F-test for the HCCM model
## Wald test
## 
## Model 1: Fv.Fm ~ 1
## Model 2: Fv.Fm ~ Topography
##   Res.Df Df      F    Pr(>F)    
## 1     23                        
## 2     22  1 26.878 3.375e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Step 4 - Multiple comparisons

• not needed for 2 levels

Step 5 - Final Diagnostic Plots not needed with Robust Method (??)

• Or you don’t need to- it’s a robust method so you don’t have to check fit??

Step 6 - Effect sizes

eta-sq = R^2 Should I be using a HCCM-corrected model below?? * Must have sums of squares to calculate effect sizes! Not included in the HCCM model output…

pander(anova_stats(Mod1, digits=3), caption = "Effect sizes - Site 1 ~ Topography")

Effect sizes - Site 1 ~ Topography (continued below)

	term	df	sumsq	meansq	statistic	p.value
Topography	Topography	1	0.143	0.143	27.99	0
Shade.Time.C	Shade.Time.C	1	0	0	0.002	0.966
…3	Residuals	21	0.107	0.005	NA	NA

Table continues below

	etasq	partial.etasq	omegasq	partial.omegasq
Topography	0.571	0.571	0.54	0.529
Shade.Time.C	0	0	-0.02	-0.043
…3	NA	NA	NA	NA

	epsilonsq	cohens.f	power
Topography	0.551	1.155	1
Shade.Time.C	-0.02	0.009	0.05
…3	NA	NA	NA

SITE 2 MODEL

Variance Test

• Heteroskedastic Topography factor
• Assumptions of regression look reasonable

D2<-subset(D, Site=="2")
boxplot((Fv.Fm) ~ Topography, data = D2)

hist(D2$Shade.Time.C)

tidy(F<-fligner.test((Fv.Fm) ~ Topography, data= D2))# unequal P = 0.0434
## # A tibble: 1 × 4
##   statistic p.value parameter method                                          
##       <dbl>   <dbl>     <dbl> <chr>                                           
## 1      6.27  0.0434         2 Fligner-Killeen test of homogeneity of variances
bptest(Fv.Fm ~ Topography*Shade.Time.C, data= D2)#not heteroskedastic w/ full model considered!
## 
##  studentized Breusch-Pagan test
## 
## data:  Fv.Fm ~ Topography * Shade.Time.C
## BP = 3.0968, df = 5, p-value = 0.6851

Fit ANOVA model & plot residual diagnostics

Fit Type I model to look at fit (normality & variance) before deciding to use HCCM! When to decide if interaction is to be included if the model isn’t a good fit, how can you decide???

Mod.2<-lm(Fv.Fm ~ Topography*Shade.Time.C, data=D2)
summary.aov(Mod.2)#ANOVA output - all three terms significant
##                         Df Sum Sq Mean Sq F value   Pr(>F)    
## Topography               2 0.0742  0.0371   3.772   0.0354 *  
## Shade.Time.C             1 0.4004  0.4004  40.717 6.61e-07 ***
## Topography:Shade.Time.C  2 0.0951  0.0476   4.836   0.0157 *  
## Residuals               28 0.2754  0.0098                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Check assumptions
par(mfrow=c(2,3))#set plotting for multiple plots
plot(Mod.2)

Regression type I

summary.lm(Mod.2)#Regression output - R2_site2<-0.6742
## 
## Call:
## lm(formula = Fv.Fm ~ Topography * Shade.Time.C, data = D2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.31958 -0.04232  0.01329  0.04906  0.14484 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               0.394388   0.029195  13.509 8.67e-14 ***
## TopographyF              -0.050812   0.043675  -1.163  0.25448    
## TopographyN               0.131160   0.042523   3.084  0.00455 ** 
## Shade.Time.C              0.003615   0.001915   1.888  0.06949 .  
## TopographyF:Shade.Time.C  0.009194   0.002988   3.077  0.00464 ** 
## TopographyN:Shade.Time.C  0.002788   0.002508   1.112  0.27562    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09917 on 28 degrees of freedom
## Multiple R-squared:  0.6742, Adjusted R-squared:  0.616 
## F-statistic: 11.59 on 5 and 28 DF,  p-value: 3.954e-06

Robust HCCM & Coefficients Test

estimatr::lm_robust function automatically includes robust standard errors.

#library(lmtest)
Cov2 <- hccm(Mod.2, type="hc3") #needs package 'car'
Mod2.HC3 <- coeftest(Mod.2, vcov.=Cov2)#needs lmtest package 
Mod2.HC3
## 
## t test of coefficients:
## 
##                            Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)               0.3943881  0.0256863 15.3540 3.637e-15 ***
## TopographyF              -0.0508122  0.0429254 -1.1837  0.246471    
## TopographyN               0.1311597  0.0431791  3.0376  0.005118 ** 
## Shade.Time.C              0.0036147  0.0017941  2.0148  0.053619 .  
## TopographyF:Shade.Time.C  0.0091941  0.0028998  3.1705  0.003668 ** 
## TopographyN:Shade.Time.C  0.0027884  0.0025114  1.1103  0.276317    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Robust ANOVA F-Test

• So, there is no omibus F-test for the regression?

Anova(Mod.2, type = "II", white.adjust = "hc3")
## Coefficient covariances computed by hccm()
## Analysis of Deviance Table (Type II tests)
## 
## Response: Fv.Fm
##                         Df       F    Pr(>F)    
## Topography               2 13.2147 9.090e-05 ***
## Shade.Time.C             1 38.7867 9.934e-07 ***
## Topography:Shade.Time.C  2  5.0788   0.01312 *  
## Residuals               28                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Test Interaction

• Fit reduced model & test if interaction improves fit significantly - Yes!

Mod.2red<-lm(Fv.Fm ~ Topography+Shade.Time.C, data = D2)
anova(Mod.2red, Mod.2)#Significantly improves fit!
## Analysis of Variance Table
## 
## Model 1: Fv.Fm ~ Topography + Shade.Time.C
## Model 2: Fv.Fm ~ Topography * Shade.Time.C
##   Res.Df     RSS Df Sum of Sq      F Pr(>F)  
## 1     30 0.37049                             
## 2     28 0.27537  2  0.095116 4.8358 0.0157 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Compare Slopes (from robst model)

Species Petal.Width.trend SE df lower.CL upper.CL setosa 0.9301727 0.6491360 144 -0.3528933 2.213239 versicolor 1.4263647 0.3459350 144 0.7425981 2.110131 virginica 0.6508306 0.2490791 144 0.1585071 1.143154

Compare slopes

pairs(m.lst) contrast estimate SE df t.ratio p.value setosa - versicolor -0.4961919 0.7355601 144 -0.675 0.7786 setosa - virginica 0.2793421 0.6952826 144 0.402 0.9149 versicolor - virginica 0.7755341 0.4262762 144 1.819 0.1669

#Anova(Mod.2, type = "II", white.adjust = "hc3")
# Obtain slopes
Mod.2$coefficients
##              (Intercept)              TopographyF              TopographyN 
##              0.394388134             -0.050812235              0.131159725 
##             Shade.Time.C TopographyF:Shade.Time.C TopographyN:Shade.Time.C 
##              0.003614697              0.009194087              0.002788356
#m.lst <- lstrends(m.interaction, "Species", var="Petal.Width")#can't find function

Wald Reduced Model R-squared F-test

• Fit reduced model & test if interaction improves fit significantly - No!

# Overall model F-test: 
Mod2.H0= lm(Fv.Fm ~ 1, data=D2)
Cov2 <- hccm(Mod.2, type="hc3") #needs package 'car'
waldtest(Mod2.H0, Mod.2, vcov=Cov2)#Overall F-test for the HCCM model
## Wald test
## 
## Model 1: Fv.Fm ~ 1
## Model 2: Fv.Fm ~ Topography * Shade.Time.C
##   Res.Df Df      F    Pr(>F)    
## 1     33                        
## 2     28  5 12.461 2.037e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Effect sizes

eta-sq = R^2 Can’t do without sums of squares-how to get HCCM-corrected SS’s?

pander(anova_stats(Mod.2, digits=3), caption = "Effect sizes - Site 2 ~ Topography x Shade Time")

Effect sizes - Site 2 ~ Topography x Shade Time (continued below)

	term	df	sumsq	meansq
Topography	Topography	2	0.074	0.037
Shade.Time.C	Shade.Time.C	1	0.4	0.4
Topography:Shade.Time.C	Topography:Shade.Time.C	2	0.095	0.048
…4	Residuals	28	0.275	0.01

Table continues below

	statistic	p.value	etasq	partial.etasq
Topography	3.772	0.035	0.088	0.212
Shade.Time.C	40.72	0	0.474	0.593
Topography:Shade.Time.C	4.836	0.016	0.113	0.257
…4	NA	NA	NA	NA

Table continues below

	omegasq	partial.omegasq	epsilonsq	cohens.f
Topography	0.064	0.14	0.065	0.519
Shade.Time.C	0.457	0.539	0.462	1.206
Topography:Shade.Time.C	0.088	0.184	0.089	0.588
…4	NA	NA	NA	NA

	power
Topography	0.687
Shade.Time.C	1
Topography:Shade.Time.C	0.799
…4	NA

SITE 3 MODEL

• Multicollinearity may be an issue here…coefficients & SE & P-values change when interaction term is dropped from the full model.
• Unequal variance requires HCCM.

Variance Test

• Heteroskedastic Topography factor
• Assumptions of regression look reasonable

D3<-subset(D, Site=="3")
boxplot((Fv.Fm) ~ Topography, data = D3)

hist(D3$Shade.Time.C)

tidy(F<-fligner.test((Fv.Fm) ~ Topography, data= D3))# unequal P = 0.00335
## # A tibble: 1 × 4
##   statistic p.value parameter method                                          
##       <dbl>   <dbl>     <dbl> <chr>                                           
## 1      11.4 0.00336         2 Fligner-Killeen test of homogeneity of variances

Fit ANOVA model & plot residual diagnostics

• Interaction n.s.
• Residual plots show heteroskedasticity & non-normality

Mod.3<-lm(Fv.Fm ~ Topography*Shade.Time.C, data=D3)
summary.aov(Mod.3)#ANOVA output - interaction n.s.
##                         Df Sum Sq Mean Sq F value   Pr(>F)    
## Topography               2 0.3078 0.15389   13.20 9.14e-05 ***
## Shade.Time.C             1 0.2225 0.22253   19.09 0.000155 ***
## Topography:Shade.Time.C  2 0.0166 0.00828    0.71 0.500098    
## Residuals               28 0.3263 0.01165                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Check assumptions
par(mfrow=c(2,3))#set plotting for multiple plots
plot(Mod.3)

Full Model - Try Robust HCCM & Coefficients Test

• Robust HCCM increases power to detect a difference in the S & N-zones!
• And increases significance of the Shade p-value (***).
• Estimates do not change, but SE’s do.
• But see next model run with the reduced model (no interaction term); confusing b/c p-value changes for all terms still in the model…suggesting multicollinearity is present between Topo & shade…?

#library(lmtest)
Cov3 <- hccm(Mod.3, type="hc3") #needs package 'car'
coeftest(Mod.3, vcov.=Cov3)#needs lmtest package 
## 
## t test of coefficients:
## 
##                            Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)               0.3110806  0.0474369  6.5578 4.136e-07 ***
## TopographyF              -0.0063137  0.0577842 -0.1093 0.9137730    
## TopographyN               0.2178991  0.1007746  2.1622 0.0392943 *  
## Shade.Time.C              0.0108881  0.0028303  3.8469 0.0006325 ***
## TopographyF:Shade.Time.C -0.0026800  0.0039173 -0.6841 0.4995080    
## TopographyN:Shade.Time.C -0.0070873  0.0043175 -1.6415 0.1118729    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Robust ANOVA - Full model Tests

• I think in hypothesis testing it makes the most sense to report the tables with all terms tested.
• Results: Topo n.s., Shade highly significant, & interaction n.s.
• Report in Table 2, but put final reduced model

Anova(Mod.3, type = "II", white.adjust = "hc3")
## Coefficient covariances computed by hccm()
## Analysis of Deviance Table (Type II tests)
## 
## Response: Fv.Fm
##                         Df       F    Pr(>F)    
## Topography               2  2.1150    0.1395    
## Shade.Time.C             1 22.6379 5.375e-05 ***
## Topography:Shade.Time.C  2  1.3529    0.2749    
## Residuals               28                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Additive Model - Robust HCCM & Coefficients Test

• Removing interaction term changes coefficients, SE’s, & P-values!
• This is a sign of multicollinearity…
• Run final reduced model with only significant terms & use Wald Test to compute overall F-test for the regression

# Overall model F-test: 
Mod3.H0= lm(Fv.Fm ~ 1, data=D3)
waldtest(Mod3.H0, Mod.3)#regular Type1 regression comparison
## Wald test
## 
## Model 1: Fv.Fm ~ 1
## Model 2: Fv.Fm ~ Topography * Shade.Time.C
##   Res.Df Df      F    Pr(>F)    
## 1     33                        
## 2     28  5 9.3845 2.441e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
waldtest(Mod3.H0, Mod3.HA=lm(Mod.3, vcov.=Cov1))#Overall F-test for the HCCM model
## Wald test
## 
## Model 1: Fv.Fm ~ 1
## Model 2: Fv.Fm ~ Topography * Shade.Time.C
##   Res.Df Df      F    Pr(>F)    
## 1     33                        
## 2     28  5 9.3845 2.441e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#library(lmtest)
Cov3 <- hccm(Mod.3, type="hc3") #needs package 'car'
coeftest(Mod.3, vcov.=Cov3)#needs lmtest package 
## 
## t test of coefficients:
## 
##                            Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)               0.3110806  0.0474369  6.5578 4.136e-07 ***
## TopographyF              -0.0063137  0.0577842 -0.1093 0.9137730    
## TopographyN               0.2178991  0.1007746  2.1622 0.0392943 *  
## Shade.Time.C              0.0108881  0.0028303  3.8469 0.0006325 ***
## TopographyF:Shade.Time.C -0.0026800  0.0039173 -0.6841 0.4995080    
## TopographyN:Shade.Time.C -0.0070873  0.0043175 -1.6415 0.1118729    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Wald Reduced Model R-squared F-test

• Fit reduced model & test if interaction improves fit significantly - No!

# Overall model F-test: 
Mod3.H0= lm(Fv.Fm ~ 1, data=D3)
Mod.3red<-lm(Fv.Fm ~Shade.Time.C, data = D3)
Cov3 <- hccm(Mod.3, type="hc3") #needs package 'car'
summary.lm(Mod.3)
## 
## Call:
## lm(formula = Fv.Fm ~ Topography * Shade.Time.C, data = D3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.28942 -0.03186  0.02661  0.06226  0.15387 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               0.311081   0.078989   3.938 0.000496 ***
## TopographyF              -0.006314   0.087461  -0.072 0.942965    
## TopographyN               0.217899   0.126637   1.721 0.096348 .  
## Shade.Time.C              0.010888   0.004520   2.409 0.022823 *  
## TopographyF:Shade.Time.C -0.002680   0.005028  -0.533 0.598242    
## TopographyN:Shade.Time.C -0.007087   0.006150  -1.152 0.258871    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.108 on 28 degrees of freedom
## Multiple R-squared:  0.6263, Adjusted R-squared:  0.5595 
## F-statistic: 9.384 on 5 and 28 DF,  p-value: 2.441e-05
#anova(Mod.3red, Mod.3)#Proof of concept, does not improve fit
Cov3 <- hccm(Mod.3red, type="hc3") #needs package 'car'
waldtest(Mod3.H0, Mod.3red, vcov=Cov3)#Overall F-test for the HCCM model
## Wald test
## 
## Model 1: Fv.Fm ~ 1
## Model 2: Fv.Fm ~ Shade.Time.C
##   Res.Df Df     F    Pr(>F)    
## 1     33                       
## 2     32  1 57.63 1.197e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Effect sizes

eta-sq = R^2 Can’t do without sums of squares-how to get HCCM-corrected SS’s?

pander(anova_stats(Mod.2, digits=3), caption = "Effect sizes - Site 2 ~ Topography x Shade Time")

Effect sizes - Site 2 ~ Topography x Shade Time (continued below)

	term	df	sumsq	meansq
Topography	Topography	2	0.074	0.037
Shade.Time.C	Shade.Time.C	1	0.4	0.4
Topography:Shade.Time.C	Topography:Shade.Time.C	2	0.095	0.048
…4	Residuals	28	0.275	0.01

Table continues below

	statistic	p.value	etasq	partial.etasq
Topography	3.772	0.035	0.088	0.212
Shade.Time.C	40.72	0	0.474	0.593
Topography:Shade.Time.C	4.836	0.016	0.113	0.257
…4	NA	NA	NA	NA

Table continues below

	omegasq	partial.omegasq	epsilonsq	cohens.f
Topography	0.064	0.14	0.065	0.519
Shade.Time.C	0.457	0.539	0.462	1.206
Topography:Shade.Time.C	0.088	0.184	0.089	0.588
…4	NA	NA	NA	NA

	power
Topography	0.687
Shade.Time.C	1
Topography:Shade.Time.C	0.799
…4	NA